Cohort 2 Taxonomy
Liz Mallott
1/21/2021
Overview of analysis
This analysis examines if overall community composition differs based on self-reported ethnicity and what specific taxa are differentially abundant. MetaPhlAn (v3.0.7) was used to profile the taxonomic composition of metagenomes and calculate unweighted and weighted UniFrac distances.
Setting up the environment
Load the necessary libraries.
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──## ✓ ggplot2 3.3.2 ✓ purrr 0.3.4
## ✓ tibble 3.0.4 ✓ dplyr 1.0.2
## ✓ tidyr 1.1.2 ✓ stringr 1.4.0
## ✓ readr 1.4.0 ✓ forcats 0.5.0## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()library(vegan)## Loading required package: permute## Loading required package: lattice## This is vegan 2.5-7library(ggplot2)
library(cowplot)
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recode## The following object is masked from 'package:purrr':
##
## somelibrary(multcomp)## Loading required package: mvtnorm## Loading required package: survival## Loading required package: TH.data## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select##
## Attaching package: 'TH.data'## The following object is masked from 'package:MASS':
##
## geyserlibrary(ggpubr)##
## Attaching package: 'ggpubr'## The following object is masked from 'package:cowplot':
##
## get_legendImporting data
Import metadata
metadata = read_tsv("year2_sample_mapping_metaphlan.tsv")##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## .default = col_character(),
## seq_id_simple = col_double(),
## fecal_use = col_double(),
## oral_use = col_double(),
## mock = col_double(),
## blank = col_double(),
## day_number = col_double(),
## time = col_time(format = ""),
## extraction_kit_use = col_double(),
## extraction = col_double(),
## library_plate = col_double(),
## id_num = col_double()
## )
## ℹ Use `spec()` for the full column specifications.Import species-level relative abundance table, transpose, and recalculate relative abundance to sum to 1 instead of 100.
relab = read_tsv("year2_merged_abundance_table_species.txt")##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## .default = col_double(),
## body_site = col_character()
## )
## ℹ Use `spec()` for the full column specifications.relab1 = relab %>% mutate_if(is.numeric, ~./100)
relab1_t = relab1 %>% gather(file, count, 2:164) %>%
spread(body_site, count)Import UniFrac matrices
wuni_fecal = as.dist(
read.table("year2_merged_abundance_wunifrac_fecal.tsv"))
uuni_fecal = as.dist(
read.table("year2_merged_abundance_uunifrac_fecal.tsv"))
wuni_oral = as.dist(
read.table("year2_merged_abundance_wunifrac_oral.tsv"))
uuni_oral = as.dist(
read.table("year2_merged_abundance_uunifrac_oral.tsv"))
wuni_fecalb = as.dist(
read.table("year2_merged_abundance_wunifrac_fecalb.tsv"))
uuni_fecalb = as.dist(
read.table("year2_merged_abundance_uunifrac_fecalb.tsv"))
wuni_oralb = as.dist(
read.table("year2_merged_abundance_wunifrac_oralb.tsv"))
uuni_oralb = as.dist(
read.table("year2_merged_abundance_uunifrac_oralb.tsv"))
wuni_fecala = as.dist(
read.table("year2_merged_abundance_wunifrac_fecala.tsv"))
uuni_fecala = as.dist(
read.table("year2_merged_abundance_uunifrac_fecala.tsv"))
wuni_orala = as.dist(
read.table("year2_merged_abundance_wunifrac_orala.tsv"))
uuni_orala = as.dist(
read.table("year2_merged_abundance_uunifrac_orala.tsv"))
wuni_fecald = as.dist(
read.table("year2_merged_abundance_wunifrac_fecald.tsv"))
uuni_fecald = as.dist(
read.table("year2_merged_abundance_uunifrac_fecald.tsv"))Filter samples
Filter out samples that have no identified species or are controls (blanks and mock communities). Separate oral samples from fecal samples.
relab_fecal = relab1_t %>% left_join(metadata, by = "file") %>%
filter(source == "fecal")
relab_fecal_filtered = relab_fecal %>%
mutate(total = rowSums(relab_fecal[,2:496])) %>%
filter(total > 0)
relab_oral = relab1_t %>% left_join(metadata, by = "file") %>%
filter(source == "oral")Calculating beta diversity metrics
UniFrac distances are calculated by a MetaPhlAn helper function. Vegan (2.5-7) is used to calculate the Bray-Curtis dissimilarity and Jaccard similarity metrics.
brayf = vegdist(relab_fecal_filtered[,2:496], method = "bray")
jaccf = vegdist(relab_fecal_filtered[,2:496], method = "jaccard")
brayo = vegdist(relab_oral[,2:496], method = "bray")
jacco = vegdist(relab_oral[,2:496], method = "jaccard")PERMANOVAs - All Samples
Differences in the taxonomic structure of the community associated with self-reported ethnicity is assessed using permutational multivariate analysis of variance.
Bray-Curtis
Fecal
adonis2(brayf ~ ethnicity, data = relab_fecal_filtered,
permutations = 4999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 4999
##
## adonis2(formula = brayf ~ ethnicity, data = relab_fecal_filtered, permutations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## ethnicity 1 1.4525 0.06536 7.6225 2e-04 ***
## Residual 109 20.7702 0.93464
## Total 110 22.2227 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Oral
adonis2(brayo ~ ethnicity, data = relab_oral,
permutations = 4999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 4999
##
## adonis2(formula = brayo ~ ethnicity, data = relab_oral, permutations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## ethnicity 1 0.3486 0.0908 3.1957 0.0042 **
## Residual 32 3.4905 0.9092
## Total 33 3.8391 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1NMDS visualization of Bray-Curtis dissimilarity with all samples.
## Run 0 stress 0.1730196
## Run 1 stress 0.2313585
## Run 2 stress 0.1900919
## Run 3 stress 0.2471194
## Run 4 stress 0.1900924
## Run 5 stress 0.2006392
## Run 6 stress 0.206551
## Run 7 stress 0.1972424
## Run 8 stress 0.2333067
## Run 9 stress 0.2111622
## Run 10 stress 0.1857487
## Run 11 stress 0.2300488
## Run 12 stress 0.173776
## Run 13 stress 0.2264283
## Run 14 stress 0.173776
## Run 15 stress 0.192742
## Run 16 stress 0.1744768
## Run 17 stress 0.2105569
## Run 18 stress 0.2032629
## Run 19 stress 0.2129299
## Run 20 stress 0.2036804
## Run 21 stress 0.2029093
## Run 22 stress 0.2086552
## Run 23 stress 0.1867749
## Run 24 stress 0.2201923
## Run 25 stress 0.1777439
## Run 26 stress 0.1730196
## ... New best solution
## ... Procrustes: rmse 1.773753e-05 max resid 8.728979e-05
## ... Similar to previous best
## *** Solution reached
## Run 0 stress 0.1592337
## Run 1 stress 0.158507
## ... New best solution
## ... Procrustes: rmse 0.009437815 max resid 0.04843542
## Run 2 stress 0.158507
## ... Procrustes: rmse 7.10493e-06 max resid 2.3685e-05
## ... Similar to previous best
## Run 3 stress 0.1869209
## Run 4 stress 0.3863172
## Run 5 stress 0.1591203
## Run 6 stress 0.1817134
## Run 7 stress 0.1660616
## Run 8 stress 0.1664842
## Run 9 stress 0.1797591
## Run 10 stress 0.1632234
## Run 11 stress 0.1632234
## Run 12 stress 0.1665592
## Run 13 stress 0.1629265
## Run 14 stress 0.1609917
## Run 15 stress 0.1798961
## Run 16 stress 0.1665592
## Run 17 stress 0.2011488
## Run 18 stress 0.1665592
## Run 19 stress 0.1632234
## Run 20 stress 0.179896
## *** Solution reached
tiff(file = "figures/nmds_plot_taxa_bray.tif", res = 150, width = 18,
height = 7, units = "in")
legend1 = get_legend(brayf + theme(legend.box.margin = margin(0, 0, 0, 12)))
col1 = plot_grid(brayf + theme(legend.position = "none"),
nrow = 1, ncol = 1, labels = c('A'),
label_size = 20)
col2 = plot_grid(brayo + theme(legend.position = "none"),
nrow = 1, ncol = 1, labels = c('B'),
label_size = 20)
plot = plot_grid(col1, col2, legend1, nrow = 1, ncol = 3,
rel_widths = c(1.5, 1.5, 0.75), align = "h", axis = "t")
plot
dev.off()## quartz_off_screen
## 2Jaccard
Fecal
adonis2(jaccf ~ ethnicity, data = relab_fecal_filtered,
permutations = 4999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 4999
##
## adonis2(formula = jaccf ~ ethnicity, data = relab_fecal_filtered, permutations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## ethnicity 1 1.904 0.05969 6.9197 2e-04 ***
## Residual 109 29.992 0.94031
## Total 110 31.896 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Oral
adonis2(jacco ~ ethnicity, data = relab_oral,
permutations = 4999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 4999
##
## adonis2(formula = jacco ~ ethnicity, data = relab_oral, permutations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## ethnicity 1 0.4911 0.07269 2.5084 0.0038 **
## Residual 32 6.2653 0.92731
## Total 33 6.7564 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1NMDS visualization of Jaccard similarity with all samples.
## Run 0 stress 0.173776
## Run 1 stress 0.1730198
## ... New best solution
## ... Procrustes: rmse 0.03083177 max resid 0.1163365
## Run 2 stress 0.2027132
## Run 3 stress 0.1977463
## Run 4 stress 0.2056424
## Run 5 stress 0.2249779
## Run 6 stress 0.2209721
## Run 7 stress 0.2121922
## Run 8 stress 0.2191132
## Run 9 stress 0.2062637
## Run 10 stress 0.1878283
## Run 11 stress 0.2035525
## Run 12 stress 0.2064486
## Run 13 stress 0.2054211
## Run 14 stress 0.2063329
## Run 15 stress 0.1777443
## Run 16 stress 0.2231832
## Run 17 stress 0.2035501
## Run 18 stress 0.2056861
## Run 19 stress 0.1884853
## Run 20 stress 0.1919388
## Run 21 stress 0.1999187
## Run 22 stress 0.2057736
## Run 23 stress 0.173776
## Run 24 stress 0.2056837
## Run 25 stress 0.1973273
## Run 26 stress 0.2022282
## Run 27 stress 0.2091878
## Run 28 stress 0.2174981
## Run 29 stress 0.2326594
## Run 30 stress 0.2082968
## Run 31 stress 0.1974395
## Run 32 stress 0.1950047
## Run 33 stress 0.1827343
## Run 34 stress 0.1930404
## Run 35 stress 0.2053464
## Run 36 stress 0.1880971
## Run 37 stress 0.1857488
## Run 38 stress 0.2166193
## Run 39 stress 0.2192174
## Run 40 stress 0.2357384
## Run 41 stress 0.2166677
## Run 42 stress 0.2171435
## Run 43 stress 0.2279408
## Run 44 stress 0.1730199
## ... Procrustes: rmse 7.797651e-05 max resid 0.0003221037
## ... Similar to previous best
## *** Solution reached
## Run 0 stress 0.1592336
## Run 1 stress 0.1624179
## Run 2 stress 0.160991
## Run 3 stress 0.2056952
## Run 4 stress 0.2032465
## Run 5 stress 0.1835446
## Run 6 stress 0.1609925
## Run 7 stress 0.2177963
## Run 8 stress 0.1930489
## Run 9 stress 0.1629265
## Run 10 stress 0.1665592
## Run 11 stress 0.2176376
## Run 12 stress 0.1622289
## Run 13 stress 0.1968443
## Run 14 stress 0.2073231
## Run 15 stress 0.1823253
## Run 16 stress 0.1676377
## Run 17 stress 0.1629265
## Run 18 stress 0.162229
## Run 19 stress 0.175451
## Run 20 stress 0.1754509
## Run 21 stress 0.1676378
## Run 22 stress 0.1591202
## ... New best solution
## ... Procrustes: rmse 0.02981526 max resid 0.1456574
## Run 23 stress 0.2009957
## Run 24 stress 0.1609921
## Run 25 stress 0.160991
## Run 26 stress 0.1859658
## Run 27 stress 0.1766694
## Run 28 stress 0.1798963
## Run 29 stress 0.1862181
## Run 30 stress 0.1629265
## Run 31 stress 0.1665592
## Run 32 stress 0.1632234
## Run 33 stress 0.2038537
## Run 34 stress 0.162229
## Run 35 stress 0.1632234
## Run 36 stress 0.1835444
## Run 37 stress 0.1609911
## Run 38 stress 0.1629265
## Run 39 stress 0.1622289
## Run 40 stress 0.1609928
## Run 41 stress 0.1699121
## Run 42 stress 0.158507
## ... New best solution
## ... Procrustes: rmse 0.03767444 max resid 0.1911456
## Run 43 stress 0.2017916
## Run 44 stress 0.2106544
## Run 45 stress 0.158507
## ... New best solution
## ... Procrustes: rmse 1.27024e-05 max resid 4.391212e-05
## ... Similar to previous best
## *** Solution reached
tiff(file = "figures/nmds_plot_taxa_jacc.tif", res = 150, width = 18,
height = 7, units = "in")
legend1 = get_legend(jaccf + theme(legend.box.margin = margin(0, 0, 0, 12)))
col1 = plot_grid(jaccf + theme(legend.position = "none"),
nrow = 1, ncol = 1, labels = c('A'),
label_size = 20)
col2 = plot_grid(jacco + theme(legend.position = "none"),
nrow = 1, ncol = 1, labels = c('B'),
label_size = 20)
plot = plot_grid(col1, col2, legend1, nrow = 1, ncol = 3,
rel_widths = c(1.5, 1.5, 0.75), align = "h", axis = "t")
plot
dev.off()## quartz_off_screen
## 2Weighted UniFrac
Fecal
adonis2(wuni_fecal ~ ethnicity, data = relab_fecal_filtered,
permutations = 4999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 4999
##
## adonis2(formula = wuni_fecal ~ ethnicity, data = relab_fecal_filtered, permutations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## ethnicity 1 0.2429 0.04942 5.6666 0.0012 **
## Residual 109 4.6721 0.95058
## Total 110 4.9150 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Oral
adonis2(wuni_oral ~ ethnicity, data = relab_oral,
permutations = 4999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 4999
##
## adonis2(formula = wuni_oral ~ ethnicity, data = relab_oral, permutations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## ethnicity 1 0.07415 0.04607 1.5456 0.1934
## Residual 32 1.53533 0.95393
## Total 33 1.60948 1.00000NMDS visualization of weighted Unifrac distances with all samples.
## Run 0 stress 0.1128885
## Run 1 stress 0.1210185
## Run 2 stress 0.1144106
## Run 3 stress 0.1144077
## Run 4 stress 0.130739
## Run 5 stress 0.1129022
## ... Procrustes: rmse 0.0009940589 max resid 0.006519917
## ... Similar to previous best
## Run 6 stress 0.1168925
## Run 7 stress 0.1128926
## ... Procrustes: rmse 0.006453419 max resid 0.0635843
## Run 8 stress 0.1128787
## ... New best solution
## ... Procrustes: rmse 0.001119216 max resid 0.009342401
## ... Similar to previous best
## Run 9 stress 0.1163883
## Run 10 stress 0.1380249
## Run 11 stress 0.1486133
## Run 12 stress 0.1425132
## Run 13 stress 0.1129707
## ... Procrustes: rmse 0.006421709 max resid 0.06307331
## Run 14 stress 0.1352977
## Run 15 stress 0.113063
## ... Procrustes: rmse 0.009634561 max resid 0.07258814
## Run 16 stress 0.1132554
## ... Procrustes: rmse 0.00704866 max resid 0.07112905
## Run 17 stress 0.1144078
## Run 18 stress 0.1129007
## ... Procrustes: rmse 0.001398501 max resid 0.009298345
## ... Similar to previous best
## Run 19 stress 0.1204263
## Run 20 stress 0.1128765
## ... New best solution
## ... Procrustes: rmse 0.0007308774 max resid 0.005945148
## ... Similar to previous best
## *** Solution reached
## Run 0 stress 0.08960119
## Run 1 stress 0.09882517
## Run 2 stress 0.1031962
## Run 3 stress 0.1126862
## Run 4 stress 0.115512
## Run 5 stress 0.1086049
## Run 6 stress 0.1073695
## Run 7 stress 0.1027136
## Run 8 stress 0.1179109
## Run 9 stress 0.1077257
## Run 10 stress 0.0896885
## ... Procrustes: rmse 0.01592665 max resid 0.07752264
## Run 11 stress 0.08960119
## ... Procrustes: rmse 2.834432e-05 max resid 0.0001151919
## ... Similar to previous best
## Run 12 stress 0.0896885
## ... Procrustes: rmse 0.01590753 max resid 0.07742581
## Run 13 stress 0.1027136
## Run 14 stress 0.1030229
## Run 15 stress 0.09882528
## Run 16 stress 0.1025333
## Run 17 stress 0.1086048
## Run 18 stress 0.08968854
## ... Procrustes: rmse 0.01593447 max resid 0.0774934
## Run 19 stress 0.1004283
## Run 20 stress 0.1085053
## *** Solution reached
tiff(file = "figures/nmds_plot_taxa_wuni.tif", res = 150, width = 18,
height = 7, units = "in")
legend1 = get_legend(wunif + theme(legend.box.margin = margin(0, 0, 0, 12)))
col1 = plot_grid(wunif + theme(legend.position = "none"),
nrow = 1, ncol = 1, labels = c('A'),
label_size = 20)
col2 = plot_grid(wunio + theme(legend.position = "none"),
nrow = 1, ncol = 1, labels = c('B'),
label_size = 20)
plot = plot_grid(col1, col2, legend1, nrow = 1, ncol = 3,
rel_widths = c(1.5, 1.5, 0.75), align = "h", axis = "t")
plot
dev.off()## quartz_off_screen
## 2Unweighted UniFrac
Fecal
adonis2(uuni_fecal ~ ethnicity, data = relab_fecal_filtered,
permutations = 4999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 4999
##
## adonis2(formula = uuni_fecal ~ ethnicity, data = relab_fecal_filtered, permutations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## ethnicity 1 0.7732 0.07238 8.5048 2e-04 ***
## Residual 109 9.9100 0.92762
## Total 110 10.6832 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Oral
adonis2(uuni_oral ~ ethnicity, data = relab_oral,
permutations = 4999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 4999
##
## adonis2(formula = uuni_oral ~ ethnicity, data = relab_oral, permutations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## ethnicity 1 0.05327 0.0402 1.3401 0.1886
## Residual 32 1.27193 0.9598
## Total 33 1.32520 1.0000NMDS visualization of Bray-Curtis dissimilarity with all samples.
## Run 0 stress 7.26794e-05
## Run 1 stress 9.69908e-05
## ... Procrustes: rmse 3.549433e-05 max resid 0.0001432463
## ... Similar to previous best
## Run 2 stress 5.901279e-05
## ... New best solution
## ... Procrustes: rmse 3.197825e-05 max resid 9.467394e-05
## ... Similar to previous best
## Run 3 stress 9.395146e-05
## ... Procrustes: rmse 4.761152e-05 max resid 0.0001204205
## ... Similar to previous best
## Run 4 stress 9.54084e-05
## ... Procrustes: rmse 2.881724e-05 max resid 8.458948e-05
## ... Similar to previous best
## Run 5 stress 7.52498e-05
## ... Procrustes: rmse 2.551549e-05 max resid 7.243461e-05
## ... Similar to previous best
## Run 6 stress 7.467214e-05
## ... Procrustes: rmse 2.395832e-05 max resid 5.006309e-05
## ... Similar to previous best
## Run 7 stress 9.933767e-05
## ... Procrustes: rmse 3.384435e-05 max resid 0.0001055956
## ... Similar to previous best
## Run 8 stress 6.422304e-05
## ... Procrustes: rmse 2.370323e-05 max resid 6.051647e-05
## ... Similar to previous best
## Run 9 stress 6.900282e-05
## ... Procrustes: rmse 2.039856e-05 max resid 4.60642e-05
## ... Similar to previous best
## Run 10 stress 9.428947e-05
## ... Procrustes: rmse 2.37946e-05 max resid 6.96226e-05
## ... Similar to previous best
## Run 11 stress 9.90301e-05
## ... Procrustes: rmse 2.512668e-05 max resid 7.256941e-05
## ... Similar to previous best
## Run 12 stress 8.85085e-05
## ... Procrustes: rmse 2.767042e-05 max resid 7.239274e-05
## ... Similar to previous best
## Run 13 stress 7.884904e-05
## ... Procrustes: rmse 1.963874e-05 max resid 6.988197e-05
## ... Similar to previous best
## Run 14 stress 9.487501e-05
## ... Procrustes: rmse 2.634122e-05 max resid 7.477966e-05
## ... Similar to previous best
## Run 15 stress 9.233808e-05
## ... Procrustes: rmse 3.330246e-05 max resid 9.539123e-05
## ... Similar to previous best
## Run 16 stress 9.531369e-05
## ... Procrustes: rmse 3.363245e-05 max resid 0.000101479
## ... Similar to previous best
## Run 17 stress 9.688617e-05
## ... Procrustes: rmse 2.839829e-05 max resid 6.846969e-05
## ... Similar to previous best
## Run 18 stress 9.723656e-05
## ... Procrustes: rmse 4.005544e-05 max resid 8.558437e-05
## ... Similar to previous best
## Run 19 stress 8.090206e-05
## ... Procrustes: rmse 2.143516e-05 max resid 6.771863e-05
## ... Similar to previous best
## Run 20 stress 8.926949e-05
## ... Procrustes: rmse 3.082847e-05 max resid 9.812705e-05
## ... Similar to previous best
## *** Solution reached## Warning in metaMDS(uuni_fecal, k = 2, trymax = 500): stress is (nearly) zero:
## you may have insufficient data## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
## Run 0 stress 0.1791554
## Run 1 stress 0.1822924
## Run 2 stress 0.1954749
## Run 3 stress 0.1806555
## Run 4 stress 0.1839397
## Run 5 stress 0.1811046
## Run 6 stress 0.1817662
## Run 7 stress 0.2109138
## Run 8 stress 0.1829339
## Run 9 stress 0.1791029
## ... New best solution
## ... Procrustes: rmse 0.007979483 max resid 0.03643087
## Run 10 stress 0.1839398
## Run 11 stress 0.1792178
## ... Procrustes: rmse 0.0194101 max resid 0.05380965
## Run 12 stress 0.1822924
## Run 13 stress 0.1825586
## Run 14 stress 0.1933016
## Run 15 stress 0.1791553
## ... Procrustes: rmse 0.007764313 max resid 0.03549337
## Run 16 stress 0.1809487
## Run 17 stress 0.186249
## Run 18 stress 0.1912585
## Run 19 stress 0.1954286
## Run 20 stress 0.1833687
## Run 21 stress 0.1829339
## Run 22 stress 0.1826122
## Run 23 stress 0.1810658
## Run 24 stress 0.1792178
## ... Procrustes: rmse 0.0193816 max resid 0.05387272
## Run 25 stress 0.1814989
## Run 26 stress 0.1819031
## Run 27 stress 0.1933562
## Run 28 stress 0.182295
## Run 29 stress 0.1791553
## ... Procrustes: rmse 0.007764574 max resid 0.03552299
## Run 30 stress 0.1787822
## ... New best solution
## ... Procrustes: rmse 0.01530204 max resid 0.05363203
## Run 31 stress 0.2137706
## Run 32 stress 0.2208541
## Run 33 stress 0.1787829
## ... Procrustes: rmse 0.0007134825 max resid 0.002664038
## ... Similar to previous best
## *** Solution reached
tiff(file = "figures/nmds_plot_taxa_uuni.tif", res = 150, width = 18,
height = 7, units = "in")
legend1 = get_legend(uunif + theme(legend.box.margin = margin(0, 0, 0, 12)))## Warning in MASS::cov.trob(data[, vars]): Probable convergence failurecol1 = plot_grid(uunif + theme(legend.position = "none"),
nrow = 1, ncol = 1, labels = c('A'),
label_size = 20)## Warning in MASS::cov.trob(data[, vars]): Probable convergence failurecol2 = plot_grid(uunio + theme(legend.position = "none"),
nrow = 1, ncol = 1, labels = c('B'),
label_size = 20)
plot = plot_grid(col1, col2, legend1, nrow = 1, ncol = 3,
rel_widths = c(1.5, 1.5, 0.75), align = "h", axis = "t")
plot
dev.off()## quartz_off_screen
## 2PERMANOVAs - Before samples
Filtering samples
Filter samples from before the diet intervention and recalculate Bray-Curtis and Jaccard metrics
relab_fecal_filtered_before = relab_fecal_filtered %>%
filter(Study_period == "Before")
relab_oral_before = relab_oral %>%
filter(Study_period == "Before")
brayfb = vegdist(relab_fecal_filtered_before[,2:496], method = "bray")
jaccfb = vegdist(relab_fecal_filtered_before[,2:496], method = "jaccard")
brayob = vegdist(relab_oral_before[,2:496], method = "bray")
jaccob = vegdist(relab_oral_before[,2:496], method = "jaccard")Bray-Curtis
Fecal
adonis2(brayfb ~ ethnicity, data = relab_fecal_filtered_before,
permutations = 4999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 4999
##
## adonis2(formula = brayfb ~ ethnicity, data = relab_fecal_filtered_before, permutations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## ethnicity 1 0.2744 0.06875 1.255 0.221
## Residual 17 3.7171 0.93125
## Total 18 3.9915 1.00000Oral
adonis2(brayob ~ ethnicity, data = relab_oral_before,
permutations = 4999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 4999
##
## adonis2(formula = brayob ~ ethnicity, data = relab_oral_before, permutations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## ethnicity 1 0.18873 0.09118 1.505 0.1206
## Residual 15 1.88112 0.90882
## Total 16 2.06985 1.00000NMDS visualization of Bray-Curtis dissimilarity with before samples.
## Run 0 stress 0.1284112
## Run 1 stress 0.1689742
## Run 2 stress 0.1284112
## ... Procrustes: rmse 1.267885e-06 max resid 2.283503e-06
## ... Similar to previous best
## Run 3 stress 0.1871944
## Run 4 stress 0.1697207
## Run 5 stress 0.1906031
## Run 6 stress 0.2872607
## Run 7 stress 0.1525329
## Run 8 stress 0.1284112
## ... Procrustes: rmse 3.129066e-06 max resid 6.608683e-06
## ... Similar to previous best
## Run 9 stress 0.1792602
## Run 10 stress 0.1522154
## Run 11 stress 0.1689742
## Run 12 stress 0.1284112
## ... Procrustes: rmse 1.700997e-05 max resid 3.624782e-05
## ... Similar to previous best
## Run 13 stress 0.1525329
## Run 14 stress 0.1792602
## Run 15 stress 0.1444672
## Run 16 stress 0.2020777
## Run 17 stress 0.1444672
## Run 18 stress 0.1284112
## ... New best solution
## ... Procrustes: rmse 1.331546e-06 max resid 2.666263e-06
## ... Similar to previous best
## Run 19 stress 0.1525331
## Run 20 stress 0.1506507
## *** Solution reached
## Run 0 stress 0.1253357
## Run 1 stress 0.1300761
## Run 2 stress 0.1384922
## Run 3 stress 0.1364623
## Run 4 stress 0.1264877
## Run 5 stress 0.1467916
## Run 6 stress 0.1294473
## Run 7 stress 0.1300761
## Run 8 stress 0.1453194
## Run 9 stress 0.1453195
## Run 10 stress 0.1264877
## Run 11 stress 0.1461133
## Run 12 stress 0.1256167
## ... Procrustes: rmse 0.02122912 max resid 0.0615593
## Run 13 stress 0.1300761
## Run 14 stress 0.1295583
## Run 15 stress 0.1253358
## ... Procrustes: rmse 0.0002664434 max resid 0.000872116
## ... Similar to previous best
## Run 16 stress 0.1455176
## Run 17 stress 0.1264877
## Run 18 stress 0.1384926
## Run 19 stress 0.1295583
## Run 20 stress 0.1364623
## *** Solution reached
tiff(file = "figures/nmds_plot_taxa_bray_before.tif", res = 150, width = 18,
height = 7, units = "in")
legend1 = get_legend(brayfb + theme(legend.box.margin = margin(0, 0, 0, 12)))
col1 = plot_grid(brayfb + theme(legend.position = "none"),
nrow = 1, ncol = 1, labels = c('A'),
label_size = 20)
col2 = plot_grid(brayob + theme(legend.position = "none"),
nrow = 1, ncol = 1, labels = c('B'),
label_size = 20)
plot = plot_grid(col1, col2, legend1, nrow = 1, ncol = 3,
rel_widths = c(1.5, 1.5, 0.75), align = "h", axis = "t")
plot
dev.off()## quartz_off_screen
## 2Jaccard
Fecal
adonis2(jaccfb ~ ethnicity, data = relab_fecal_filtered_before,
permutations = 4999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 4999
##
## adonis2(formula = jaccfb ~ ethnicity, data = relab_fecal_filtered_before, permutations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## ethnicity 1 0.4122 0.07421 1.3628 0.1418
## Residual 17 5.1416 0.92579
## Total 18 5.5537 1.00000Oral
adonis2(jaccob ~ ethnicity, data = relab_oral_before,
permutations = 4999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 4999
##
## adonis2(formula = jaccob ~ ethnicity, data = relab_oral_before, permutations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## ethnicity 1 0.2815 0.07968 1.2988 0.1374
## Residual 15 3.2508 0.92032
## Total 16 3.5323 1.00000NMDS visualization of Jaccard similarity with before samples.
## Run 0 stress 0.1284112
## Run 1 stress 0.1284112
## ... Procrustes: rmse 5.026012e-06 max resid 1.072812e-05
## ... Similar to previous best
## Run 2 stress 0.1284112
## ... Procrustes: rmse 2.255768e-06 max resid 4.757595e-06
## ... Similar to previous best
## Run 3 stress 0.1643886
## Run 4 stress 0.1689742
## Run 5 stress 0.1444672
## Run 6 stress 0.1284112
## ... Procrustes: rmse 3.127735e-06 max resid 6.192814e-06
## ... Similar to previous best
## Run 7 stress 0.1525329
## Run 8 stress 0.1284112
## ... Procrustes: rmse 8.96628e-07 max resid 1.623914e-06
## ... Similar to previous best
## Run 9 stress 0.194464
## Run 10 stress 0.152533
## Run 11 stress 0.1697207
## Run 12 stress 0.1284112
## ... New best solution
## ... Procrustes: rmse 4.819807e-07 max resid 9.644547e-07
## ... Similar to previous best
## Run 13 stress 0.1610714
## Run 14 stress 0.1842235
## Run 15 stress 0.2633195
## Run 16 stress 0.1842236
## Run 17 stress 0.1444672
## Run 18 stress 0.152533
## Run 19 stress 0.1689742
## Run 20 stress 0.1945795
## *** Solution reached
## Run 0 stress 0.1253358
## Run 1 stress 0.1481009
## Run 2 stress 0.1498922
## Run 3 stress 0.1585965
## Run 4 stress 0.1453194
## Run 5 stress 0.1463235
## Run 6 stress 0.1384927
## Run 7 stress 0.1455176
## Run 8 stress 0.1264877
## Run 9 stress 0.1463234
## Run 10 stress 0.1461133
## Run 11 stress 0.1467916
## Run 12 stress 0.1455177
## Run 13 stress 0.1453195
## Run 14 stress 0.1467915
## Run 15 stress 0.1253357
## ... New best solution
## ... Procrustes: rmse 0.0001619306 max resid 0.0005307169
## ... Similar to previous best
## Run 16 stress 0.1453195
## Run 17 stress 0.1253357
## ... Procrustes: rmse 0.000132518 max resid 0.000433773
## ... Similar to previous best
## Run 18 stress 0.1467916
## Run 19 stress 0.1461134
## Run 20 stress 0.1453194
## *** Solution reached
tiff(file = "figures/nmds_plot_taxa_jacc_before.tif", res = 150, width = 18,
height = 7, units = "in")
legend1 = get_legend(jaccfb + theme(legend.box.margin = margin(0, 0, 0, 12)))
col1 = plot_grid(jaccfb + theme(legend.position = "none"),
nrow = 1, ncol = 1, labels = c('A'),
label_size = 20)
col2 = plot_grid(jaccob + theme(legend.position = "none"),
nrow = 1, ncol = 1, labels = c('B'),
label_size = 20)
plot = plot_grid(col1, col2, legend1, nrow = 1, ncol = 3,
rel_widths = c(1.5, 1.5, 0.75), align = "h", axis = "t")
plot
dev.off()## quartz_off_screen
## 2Weighted UniFrac
Fecal
adonis2(wuni_fecalb ~ ethnicity, data = relab_fecal_filtered_before,
permutations = 4999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 4999
##
## adonis2(formula = wuni_fecalb ~ ethnicity, data = relab_fecal_filtered_before, permutations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## ethnicity 1 0.02942 0.03414 0.6009 0.6462
## Residual 17 0.83220 0.96586
## Total 18 0.86162 1.00000Oral
adonis2(wuni_oralb ~ ethnicity, data = relab_oral_before,
permutations = 4999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 4999
##
## adonis2(formula = wuni_oralb ~ ethnicity, data = relab_oral_before, permutations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## ethnicity 1 0.04241 0.04763 0.7501 0.4672
## Residual 15 0.84808 0.95237
## Total 16 0.89049 1.00000NMDS visualization of weighted Unifrac distances with before samples.
## Run 0 stress 0.07991057
## Run 1 stress 0.08273946
## Run 2 stress 0.08273945
## Run 3 stress 0.07999749
## ... Procrustes: rmse 0.009259724 max resid 0.02902095
## Run 4 stress 0.08403418
## Run 5 stress 0.08148165
## Run 6 stress 0.08148163
## Run 7 stress 0.08403424
## Run 8 stress 0.08273944
## Run 9 stress 0.07991038
## ... New best solution
## ... Procrustes: rmse 0.000109443 max resid 0.0003392373
## ... Similar to previous best
## Run 10 stress 0.08403433
## Run 11 stress 0.07999758
## ... Procrustes: rmse 0.009292403 max resid 0.02906161
## Run 12 stress 0.08148165
## Run 13 stress 0.2041272
## Run 14 stress 0.1479628
## Run 15 stress 0.1419221
## Run 16 stress 0.2470987
## Run 17 stress 0.2321205
## Run 18 stress 0.08403426
## Run 19 stress 0.08273946
## Run 20 stress 0.08148165
## *** Solution reached
## Run 0 stress 0.06894997
## Run 1 stress 0.06966296
## Run 2 stress 0.06894997
## ... Procrustes: rmse 9.138811e-06 max resid 2.108204e-05
## ... Similar to previous best
## Run 3 stress 0.07415722
## Run 4 stress 0.07643535
## Run 5 stress 0.07415714
## Run 6 stress 0.07616554
## Run 7 stress 0.07415714
## Run 8 stress 0.07415716
## Run 9 stress 0.07615363
## Run 10 stress 0.07415714
## Run 11 stress 0.06894999
## ... Procrustes: rmse 2.244914e-05 max resid 5.319055e-05
## ... Similar to previous best
## Run 12 stress 0.07415714
## Run 13 stress 0.06894997
## ... Procrustes: rmse 5.549988e-05 max resid 0.0001784903
## ... Similar to previous best
## Run 14 stress 0.06894997
## ... Procrustes: rmse 1.116434e-05 max resid 3.268703e-05
## ... Similar to previous best
## Run 15 stress 0.07415714
## Run 16 stress 0.06894997
## ... Procrustes: rmse 6.667681e-06 max resid 1.536305e-05
## ... Similar to previous best
## Run 17 stress 0.06966296
## Run 18 stress 0.07047077
## Run 19 stress 0.06895003
## ... Procrustes: rmse 0.0001233991 max resid 0.0003953277
## ... Similar to previous best
## Run 20 stress 0.06894997
## ... New best solution
## ... Procrustes: rmse 2.332943e-05 max resid 7.447102e-05
## ... Similar to previous best
## *** Solution reached
tiff(file = "figures/nmds_plot_taxa_wuni_before.tif", res = 150, width = 18,
height = 7, units = "in")
legend1 = get_legend(wunifb + theme(legend.box.margin = margin(0, 0, 0, 12)))
col1 = plot_grid(wunifb + theme(legend.position = "none"),
nrow = 1, ncol = 1, labels = c('A'),
label_size = 20)
col2 = plot_grid(wuniob + theme(legend.position = "none"),
nrow = 1, ncol = 1, labels = c('B'),
label_size = 20)
plot = plot_grid(col1, col2, legend1, nrow = 1, ncol = 3,
rel_widths = c(1.5, 1.5, 0.75), align = "h", axis = "t")
plot
dev.off()## quartz_off_screen
## 2Unweighted UniFrac
Fecal
adonis2(uuni_fecalb ~ ethnicity, data = relab_fecal_filtered_before,
permutations = 4999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 4999
##
## adonis2(formula = uuni_fecalb ~ ethnicity, data = relab_fecal_filtered_before, permutations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## ethnicity 1 0.18751 0.08793 1.6389 0.0192 *
## Residual 17 1.94502 0.91207
## Total 18 2.13253 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Oral
adonis2(uuni_oralb ~ ethnicity, data = relab_oral_before,
permutations = 4999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 4999
##
## adonis2(formula = uuni_oralb ~ ethnicity, data = relab_oral_before, permutations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## ethnicity 1 0.03204 0.04827 0.7608 0.7082
## Residual 15 0.63162 0.95173
## Total 16 0.66365 1.00000NMDS visualization of unweighted UniFrac distances with before samples.
## Run 0 stress 9.874988e-05
## Run 1 stress 9.738051e-05
## ... New best solution
## ... Procrustes: rmse 5.287388e-05 max resid 0.0001821857
## ... Similar to previous best
## Run 2 stress 9.567545e-05
## ... New best solution
## ... Procrustes: rmse 0.0001116539 max resid 0.0002738316
## ... Similar to previous best
## Run 3 stress 9.531762e-05
## ... New best solution
## ... Procrustes: rmse 0.0001555189 max resid 0.0003845121
## ... Similar to previous best
## Run 4 stress 9.671068e-05
## ... Procrustes: rmse 0.0001344295 max resid 0.0003474037
## ... Similar to previous best
## Run 5 stress 5.434395e-05
## ... New best solution
## ... Procrustes: rmse 0.000100357 max resid 0.0002018783
## ... Similar to previous best
## Run 6 stress 9.900469e-05
## ... Procrustes: rmse 0.0001047454 max resid 0.0001988247
## ... Similar to previous best
## Run 7 stress 9.291372e-05
## ... Procrustes: rmse 9.765583e-05 max resid 0.0002577777
## ... Similar to previous best
## Run 8 stress 8.722854e-05
## ... Procrustes: rmse 7.334147e-05 max resid 0.0001892078
## ... Similar to previous best
## Run 9 stress 7.584566e-05
## ... Procrustes: rmse 5.919929e-05 max resid 0.0001497131
## ... Similar to previous best
## Run 10 stress 9.377897e-05
## ... Procrustes: rmse 0.0001069664 max resid 0.0002615232
## ... Similar to previous best
## Run 11 stress 9.360577e-05
## ... Procrustes: rmse 0.0001182038 max resid 0.0003032233
## ... Similar to previous best
## Run 12 stress 9.431176e-05
## ... Procrustes: rmse 7.99195e-05 max resid 0.0001745545
## ... Similar to previous best
## Run 13 stress 9.616363e-05
## ... Procrustes: rmse 0.0001205418 max resid 0.0003031172
## ... Similar to previous best
## Run 14 stress 6.990469e-05
## ... Procrustes: rmse 5.261081e-05 max resid 0.0001056664
## ... Similar to previous best
## Run 15 stress 9.423776e-05
## ... Procrustes: rmse 9.508213e-05 max resid 0.0001781795
## ... Similar to previous best
## Run 16 stress 9.876871e-05
## ... Procrustes: rmse 8.801927e-05 max resid 0.0002047066
## ... Similar to previous best
## Run 17 stress 9.25114e-05
## ... Procrustes: rmse 7.214541e-05 max resid 0.0001535098
## ... Similar to previous best
## Run 18 stress 9.256526e-05
## ... Procrustes: rmse 0.0001124195 max resid 0.000271437
## ... Similar to previous best
## Run 19 stress 8.026396e-05
## ... Procrustes: rmse 7.486417e-05 max resid 0.0001550846
## ... Similar to previous best
## Run 20 stress 8.552818e-05
## ... Procrustes: rmse 9.18725e-05 max resid 0.0002441121
## ... Similar to previous best
## *** Solution reached## Warning in metaMDS(uuni_fecalb, k = 2, trymax = 500): stress is (nearly) zero:
## you may have insufficient data## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
## Run 0 stress 0.1694422
## Run 1 stress 0.1715766
## Run 2 stress 0.1703132
## Run 3 stress 0.2433225
## Run 4 stress 0.1694423
## ... Procrustes: rmse 6.586859e-05 max resid 0.0001790191
## ... Similar to previous best
## Run 5 stress 0.1715766
## Run 6 stress 0.1694422
## ... Procrustes: rmse 1.201997e-05 max resid 3.246919e-05
## ... Similar to previous best
## Run 7 stress 0.2186982
## Run 8 stress 0.1768019
## Run 9 stress 0.1715766
## Run 10 stress 0.1694423
## ... Procrustes: rmse 4.487855e-05 max resid 0.0001187343
## ... Similar to previous best
## Run 11 stress 0.2136743
## Run 12 stress 0.1768019
## Run 13 stress 0.1694423
## ... Procrustes: rmse 3.643282e-05 max resid 0.0001001615
## ... Similar to previous best
## Run 14 stress 0.1715766
## Run 15 stress 0.1696274
## ... Procrustes: rmse 0.01855263 max resid 0.05428006
## Run 16 stress 0.1694423
## ... Procrustes: rmse 0.0001232676 max resid 0.0003308309
## ... Similar to previous best
## Run 17 stress 0.1715766
## Run 18 stress 0.2199095
## Run 19 stress 0.1715766
## Run 20 stress 0.2198418
## *** Solution reached
tiff(file = "figures/nmds_plot_taxa_uuni_before.tif", res = 150, width = 18,
height = 7, units = "in")
legend1 = get_legend(uunifb + theme(legend.box.margin = margin(0, 0, 0, 12)))## Warning in MASS::cov.trob(data[, vars]): Probable convergence failurecol1 = plot_grid(uunifb + theme(legend.position = "none"),
nrow = 1, ncol = 1, labels = c('A'),
label_size = 20)## Warning in MASS::cov.trob(data[, vars]): Probable convergence failurecol2 = plot_grid(uuniob + theme(legend.position = "none"),
nrow = 1, ncol = 1, labels = c('B'),
label_size = 20)
plot = plot_grid(col1, col2, legend1, nrow = 1, ncol = 3,
rel_widths = c(1.5, 1.5, 0.75), align = "h", axis = "t")
plot
dev.off()## quartz_off_screen
## 2PERMANOVAs - During samples
Filter samples
Filter samples from during the diet intervention and recalculate Bray-Curtis and Jaccard metrics
relab_fecal_filtered_during = relab_fecal_filtered %>%
filter(Study_period == "During")
brayfd = vegdist(relab_fecal_filtered_during[,2:496], method = "bray")
jaccfd = vegdist(relab_fecal_filtered_during[,2:496], method = "jaccard")Bray-Curtis
Fecal
adonis2(brayfd ~ ethnicity, data = relab_fecal_filtered_during,
permutations = 4999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 4999
##
## adonis2(formula = brayfd ~ ethnicity, data = relab_fecal_filtered_during, permutations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## ethnicity 1 0.7899 0.06703 4.2387 4e-04 ***
## Residual 59 10.9957 0.93297
## Total 60 11.7856 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1NMDS visualization of Bray-Curtis dissimilarity with during samples.
## Run 0 stress 0.1664788
## Run 1 stress 0.1752781
## Run 2 stress 0.2053094
## Run 3 stress 0.1998895
## Run 4 stress 0.1752755
## Run 5 stress 0.1969035
## Run 6 stress 0.2004216
## Run 7 stress 0.2313365
## Run 8 stress 0.1981598
## Run 9 stress 0.2247829
## Run 10 stress 0.2112386
## Run 11 stress 0.217677
## Run 12 stress 0.224643
## Run 13 stress 0.2069444
## Run 14 stress 0.2063212
## Run 15 stress 0.1980561
## Run 16 stress 0.1958214
## Run 17 stress 0.1664617
## ... New best solution
## ... Procrustes: rmse 0.001075342 max resid 0.006206652
## ... Similar to previous best
## Run 18 stress 0.1982036
## Run 19 stress 0.2112386
## Run 20 stress 0.1724189
## *** Solution reached
tiff(file = "figures/nmds_plot_taxa_bray_during.tif", res = 150, width = 11,
height = 7, units = "in")
brayfd
dev.off()## quartz_off_screen
## 2Jaccard
Fecal
adonis2(jaccfd ~ ethnicity, data = relab_fecal_filtered_during,
permutations = 4999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 4999
##
## adonis2(formula = jaccfd ~ ethnicity, data = relab_fecal_filtered_during, permutations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## ethnicity 1 1.0741 0.06259 3.9392 2e-04 ***
## Residual 59 16.0879 0.93741
## Total 60 17.1620 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1NMDS visualization of Jaccard similarity with during samples.
## Run 0 stress 0.1664788
## Run 1 stress 0.2118375
## Run 2 stress 0.1724188
## Run 3 stress 0.175278
## Run 4 stress 0.2004768
## Run 5 stress 0.1664617
## ... New best solution
## ... Procrustes: rmse 0.00107627 max resid 0.006162014
## ... Similar to previous best
## Run 6 stress 0.1970234
## Run 7 stress 0.1664788
## ... Procrustes: rmse 0.0010852 max resid 0.006259243
## ... Similar to previous best
## Run 8 stress 0.2210808
## Run 9 stress 0.20746
## Run 10 stress 0.1806977
## Run 11 stress 0.1664617
## ... New best solution
## ... Procrustes: rmse 3.550179e-05 max resid 0.0001502291
## ... Similar to previous best
## Run 12 stress 0.2102412
## Run 13 stress 0.1664617
## ... Procrustes: rmse 1.514577e-05 max resid 4.866654e-05
## ... Similar to previous best
## Run 14 stress 0.2057964
## Run 15 stress 0.1998369
## Run 16 stress 0.2055727
## Run 17 stress 0.2008209
## Run 18 stress 0.175278
## Run 19 stress 0.2062881
## Run 20 stress 0.1945864
## *** Solution reached
tiff(file = "figures/nmds_plot_taxa_jacc_during.tif", res = 150, width = 11,
height = 7, units = "in")
jaccfd
dev.off()## quartz_off_screen
## 2Weighted UniFrac
Fecal
adonis2(wuni_fecald ~ ethnicity, data = relab_fecal_filtered_during,
permutations = 4999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 4999
##
## adonis2(formula = wuni_fecald ~ ethnicity, data = relab_fecal_filtered_during, permutations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## ethnicity 1 0.15963 0.06033 3.7881 0.009 **
## Residual 59 2.48636 0.93967
## Total 60 2.64599 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1NMDS visualization of weighted Unifrac distances with during samples.
## Run 0 stress 0.1065766
## Run 1 stress 0.1066787
## ... Procrustes: rmse 0.01085342 max resid 0.08224327
## Run 2 stress 0.1506047
## Run 3 stress 0.1159546
## Run 4 stress 0.1167693
## Run 5 stress 0.1164082
## Run 6 stress 0.1412151
## Run 7 stress 0.116769
## Run 8 stress 0.1168791
## Run 9 stress 0.1168802
## Run 10 stress 0.143699
## Run 11 stress 0.1182275
## Run 12 stress 0.1157516
## Run 13 stress 0.1363755
## Run 14 stress 0.1182267
## Run 15 stress 0.1405206
## Run 16 stress 0.1125866
## Run 17 stress 0.1103593
## Run 18 stress 0.1459305
## Run 19 stress 0.1186453
## Run 20 stress 0.1371618
## Run 21 stress 0.1103092
## Run 22 stress 0.1103094
## Run 23 stress 0.1066172
## ... Procrustes: rmse 0.004435014 max resid 0.02910857
## Run 24 stress 0.1066786
## ... Procrustes: rmse 0.01084442 max resid 0.08216484
## Run 25 stress 0.1159532
## Run 26 stress 0.1170446
## Run 27 stress 0.1404304
## Run 28 stress 0.1159534
## Run 29 stress 0.1066787
## ... Procrustes: rmse 0.0108583 max resid 0.08228771
## Run 30 stress 0.1339132
## Run 31 stress 0.1295888
## Run 32 stress 0.1170402
## Run 33 stress 0.1066172
## ... Procrustes: rmse 0.004433051 max resid 0.02912203
## Run 34 stress 0.1065764
## ... New best solution
## ... Procrustes: rmse 8.767366e-05 max resid 0.0004006176
## ... Similar to previous best
## *** Solution reached
tiff(file = "figures/nmds_plot_taxa_wuni_during.tif", res = 150, width = 11,
height = 7, units = "in")
wunif
dev.off()## quartz_off_screen
## 2Unweighted UniFrac
Fecal
adonis2(uuni_fecald ~ ethnicity, data = relab_fecal_filtered_during,
permutations = 4999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 4999
##
## adonis2(formula = uuni_fecald ~ ethnicity, data = relab_fecal_filtered_during, permutations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## ethnicity 1 0.4436 0.07779 4.9764 2e-04 ***
## Residual 59 5.2595 0.92221
## Total 60 5.7031 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1NMDS visualization of unweighted Unifrac distance with during samples.
## Run 0 stress 0.1793941
## Run 1 stress 0.2142523
## Run 2 stress 0.2151973
## Run 3 stress 0.1848087
## Run 4 stress 0.1793941
## ... New best solution
## ... Procrustes: rmse 7.947352e-05 max resid 0.0004143951
## ... Similar to previous best
## Run 5 stress 0.1847541
## Run 6 stress 0.2180861
## Run 7 stress 0.1815689
## Run 8 stress 0.1794329
## ... Procrustes: rmse 0.006333638 max resid 0.03047401
## Run 9 stress 0.1862549
## Run 10 stress 0.1800931
## Run 11 stress 0.1869278
## Run 12 stress 0.1796544
## ... Procrustes: rmse 0.01109283 max resid 0.03652654
## Run 13 stress 0.218491
## Run 14 stress 0.4059784
## Run 15 stress 0.4064148
## Run 16 stress 0.181569
## Run 17 stress 0.1847227
## Run 18 stress 0.2179301
## Run 19 stress 0.2156235
## Run 20 stress 0.1815697
## *** Solution reached
tiff(file = "figures/nmds_plot_taxa_uuni_during.tif", res = 150, width = 18,
height = 7, units = "in")
uunifd
dev.off()## quartz_off_screen
## 2PERMANOVAs - After samples
Filter samples
Filter after samples and recalculate Bray-Curtis and Jaccard metrics
relab_fecal_filtered_after = relab_fecal_filtered %>%
filter(Study_period == "After")
relab_oral_after = relab_oral %>%
filter(Study_period == "After")
brayfa = vegdist(relab_fecal_filtered_after[,2:496], method = "bray")
jaccfa = vegdist(relab_fecal_filtered_after[,2:496], method = "jaccard")
brayoa = vegdist(relab_oral_after[,2:496], method = "bray")
jaccoa = vegdist(relab_oral_after[,2:496], method = "jaccard")Bray-Curtis
Fecal
adonis2(brayfa ~ ethnicity, data = relab_fecal_filtered_after,
permutations = 4999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 4999
##
## adonis2(formula = brayfa ~ ethnicity, data = relab_fecal_filtered_after, permutations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## ethnicity 1 0.4762 0.07579 2.3782 0.0226 *
## Residual 29 5.8065 0.92421
## Total 30 6.2827 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Oral
adonis2(brayoa ~ ethnicity, data = relab_oral_after,
permutations = 4999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 4999
##
## adonis2(formula = brayoa ~ ethnicity, data = relab_oral_after, permutations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## ethnicity 1 0.24789 0.14292 2.5014 0.0174 *
## Residual 15 1.48655 0.85708
## Total 16 1.73444 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1NMDS visualization of Bray-Curtis dissimilarity with after samples.
## Run 0 stress 0.1590347
## Run 1 stress 0.1944615
## Run 2 stress 0.1937004
## Run 3 stress 0.1895439
## Run 4 stress 0.205105
## Run 5 stress 0.1658979
## Run 6 stress 0.1713295
## Run 7 stress 0.1711856
## Run 8 stress 0.1590348
## ... Procrustes: rmse 6.513386e-05 max resid 0.0001955462
## ... Similar to previous best
## Run 9 stress 0.2038524
## Run 10 stress 0.1714717
## Run 11 stress 0.1590347
## ... Procrustes: rmse 0.0001124624 max resid 0.0003000413
## ... Similar to previous best
## Run 12 stress 0.1658979
## Run 13 stress 0.1590347
## ... Procrustes: rmse 9.894018e-05 max resid 0.000277995
## ... Similar to previous best
## Run 14 stress 0.1866295
## Run 15 stress 0.2021307
## Run 16 stress 0.2031112
## Run 17 stress 0.1656963
## Run 18 stress 0.1658979
## Run 19 stress 0.1632431
## Run 20 stress 0.196653
## *** Solution reached
## Run 0 stress 0.1367903
## Run 1 stress 0.1635865
## Run 2 stress 0.1637332
## Run 3 stress 0.1367904
## ... Procrustes: rmse 7.20447e-05 max resid 0.0001608161
## ... Similar to previous best
## Run 4 stress 0.1504065
## Run 5 stress 0.1504066
## Run 6 stress 0.1367903
## ... Procrustes: rmse 3.209441e-05 max resid 8.068756e-05
## ... Similar to previous best
## Run 7 stress 0.1617068
## Run 8 stress 0.1345672
## ... New best solution
## ... Procrustes: rmse 0.0376222 max resid 0.1126213
## Run 9 stress 0.1769398
## Run 10 stress 0.1345672
## ... Procrustes: rmse 2.563029e-05 max resid 5.312082e-05
## ... Similar to previous best
## Run 11 stress 0.1367903
## Run 12 stress 0.1345672
## ... Procrustes: rmse 1.463215e-05 max resid 3.969615e-05
## ... Similar to previous best
## Run 13 stress 0.1345672
## ... Procrustes: rmse 1.981793e-05 max resid 4.80818e-05
## ... Similar to previous best
## Run 14 stress 0.1519729
## Run 15 stress 0.1617068
## Run 16 stress 0.1635865
## Run 17 stress 0.1617068
## Run 18 stress 0.1367905
## Run 19 stress 0.1563232
## Run 20 stress 0.1520582
## *** Solution reached
tiff(file = "figures/nmds_plot_taxa_bray_after.tif", res = 150, width = 18,
height = 7, units = "in")
legend1 = get_legend(brayfa + theme(legend.box.margin = margin(0, 0, 0, 12)))
col1 = plot_grid(brayfa + theme(legend.position = "none"),
nrow = 1, ncol = 1, labels = c('A'),
label_size = 20)
col2 = plot_grid(brayoa + theme(legend.position = "none"),
nrow = 1, ncol = 1, labels = c('B'),
label_size = 20)
plot = plot_grid(col1, col2, legend1, nrow = 1, ncol = 3,
rel_widths = c(1.5, 1.5, 0.75), align = "h", axis = "t")
plot
dev.off()## quartz_off_screen
## 2Jaccard
Fecal
adonis2(jaccfa ~ ethnicity, data = relab_fecal_filtered_after,
permutations = 4999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 4999
##
## adonis2(formula = jaccfa ~ ethnicity, data = relab_fecal_filtered_after, permutations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## ethnicity 1 0.6328 0.07096 2.2149 0.0156 *
## Residual 29 8.2848 0.92904
## Total 30 8.9175 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Oral
adonis2(jaccoa ~ ethnicity, data = relab_oral_after,
permutations = 4999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 4999
##
## adonis2(formula = jaccoa ~ ethnicity, data = relab_oral_after, permutations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## ethnicity 1 0.36814 0.11723 1.9919 0.0196 *
## Residual 15 2.77222 0.88277
## Total 16 3.14036 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1NMDS visualization of Jaccard similarity with after samples.
## Run 0 stress 0.1590347
## Run 1 stress 0.1658979
## Run 2 stress 0.2130876
## Run 3 stress 0.2118508
## Run 4 stress 0.2053138
## Run 5 stress 0.1871108
## Run 6 stress 0.1712887
## Run 7 stress 0.1590348
## ... Procrustes: rmse 0.0001206207 max resid 0.0003348898
## ... Similar to previous best
## Run 8 stress 0.1714715
## Run 9 stress 0.1632431
## Run 10 stress 0.1714068
## Run 11 stress 0.1926743
## Run 12 stress 0.1773822
## Run 13 stress 0.1866586
## Run 14 stress 0.1590347
## ... Procrustes: rmse 3.398221e-05 max resid 9.946736e-05
## ... Similar to previous best
## Run 15 stress 0.1816165
## Run 16 stress 0.1867049
## Run 17 stress 0.1590348
## ... Procrustes: rmse 5.010818e-05 max resid 0.0001729842
## ... Similar to previous best
## Run 18 stress 0.1656963
## Run 19 stress 0.1590348
## ... Procrustes: rmse 0.0001394517 max resid 0.0004449221
## ... Similar to previous best
## Run 20 stress 0.1912909
## *** Solution reached
## Run 0 stress 0.1367902
## Run 1 stress 0.1504065
## Run 2 stress 0.1617068
## Run 3 stress 0.1345672
## ... New best solution
## ... Procrustes: rmse 0.03748586 max resid 0.1125334
## Run 4 stress 0.1367902
## Run 5 stress 0.1563232
## Run 6 stress 0.1563232
## Run 7 stress 0.1367909
## Run 8 stress 0.1367905
## Run 9 stress 0.1637331
## Run 10 stress 0.1367902
## Run 11 stress 0.1577926
## Run 12 stress 0.1637332
## Run 13 stress 0.1367907
## Run 14 stress 0.1367904
## Run 15 stress 0.1769394
## Run 16 stress 0.1367902
## Run 17 stress 0.1367904
## Run 18 stress 0.1563232
## Run 19 stress 0.1345672
## ... Procrustes: rmse 2.443653e-05 max resid 4.620881e-05
## ... Similar to previous best
## Run 20 stress 0.1345672
## ... New best solution
## ... Procrustes: rmse 2.487242e-05 max resid 4.789656e-05
## ... Similar to previous best
## *** Solution reached
tiff(file = "figures/nmds_plot_taxa_jacc_after.tif", res = 150, width = 18,
height = 7, units = "in")
legend1 = get_legend(jaccfa + theme(legend.box.margin = margin(0, 0, 0, 12)))
col1 = plot_grid(jaccfa + theme(legend.position = "none"),
nrow = 1, ncol = 1, labels = c('A'),
label_size = 20)
col2 = plot_grid(jaccoa + theme(legend.position = "none"),
nrow = 1, ncol = 1, labels = c('B'),
label_size = 20)
plot = plot_grid(col1, col2, legend1, nrow = 1, ncol = 3,
rel_widths = c(1.5, 1.5, 0.75), align = "h", axis = "t")
plot
dev.off()## quartz_off_screen
## 2Weighted UniFrac
Fecal
adonis2(wuni_fecala ~ ethnicity, data = relab_fecal_filtered_after,
permutations = 4999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 4999
##
## adonis2(formula = wuni_fecala ~ ethnicity, data = relab_fecal_filtered_after, permutations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## ethnicity 1 0.09394 0.07012 2.1868 0.0774 .
## Residual 29 1.24577 0.92988
## Total 30 1.33971 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Oral
adonis2(wuni_orala ~ ethnicity, data = relab_oral_after,
permutations = 4999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 4999
##
## adonis2(formula = wuni_orala ~ ethnicity, data = relab_oral_after, permutations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## ethnicity 1 0.07966 0.11412 1.9323 0.1236
## Residual 15 0.61840 0.88588
## Total 16 0.69806 1.00000NMDS visualization of weighted Unifrac distances with after samples.
## Run 0 stress 0.103874
## Run 1 stress 0.1030343
## ... New best solution
## ... Procrustes: rmse 0.03049715 max resid 0.1471977
## Run 2 stress 0.127867
## Run 3 stress 0.1283418
## Run 4 stress 0.1340227
## Run 5 stress 0.1108314
## Run 6 stress 0.1079067
## Run 7 stress 0.1079067
## Run 8 stress 0.1079067
## Run 9 stress 0.103874
## Run 10 stress 0.1038767
## Run 11 stress 0.1254346
## Run 12 stress 0.1079067
## Run 13 stress 0.1079068
## Run 14 stress 0.1079068
## Run 15 stress 0.1079068
## Run 16 stress 0.112468
## Run 17 stress 0.1079067
## Run 18 stress 0.1328807
## Run 19 stress 0.103874
## Run 20 stress 0.1079068
## Run 21 stress 0.1079067
## Run 22 stress 0.1038762
## Run 23 stress 0.1030343
## ... New best solution
## ... Procrustes: rmse 8.352202e-05 max resid 0.0003853117
## ... Similar to previous best
## *** Solution reached
## Run 0 stress 0.05926746
## Run 1 stress 0.07191275
## Run 2 stress 0.08726035
## Run 3 stress 0.07191275
## Run 4 stress 0.05926746
## ... Procrustes: rmse 1.319554e-06 max resid 4.077547e-06
## ... Similar to previous best
## Run 5 stress 0.08358197
## Run 6 stress 0.05926746
## ... Procrustes: rmse 1.720742e-06 max resid 4.009741e-06
## ... Similar to previous best
## Run 7 stress 0.05926746
## ... New best solution
## ... Procrustes: rmse 4.973794e-06 max resid 1.59931e-05
## ... Similar to previous best
## Run 8 stress 0.07401827
## Run 9 stress 0.07401828
## Run 10 stress 0.09924667
## Run 11 stress 0.07401827
## Run 12 stress 0.07191275
## Run 13 stress 0.07191275
## Run 14 stress 0.07191275
## Run 15 stress 0.08605518
## Run 16 stress 0.09687146
## Run 17 stress 0.09478662
## Run 18 stress 0.07401827
## Run 19 stress 0.05926746
## ... New best solution
## ... Procrustes: rmse 9.301701e-07 max resid 2.343873e-06
## ... Similar to previous best
## Run 20 stress 0.07191275
## *** Solution reached
tiff(file = "figures/nmds_plot_taxa_wuni_after.tif", res = 150, width = 18,
height = 7, units = "in")
legend1 = get_legend(wunifa + theme(legend.box.margin = margin(0, 0, 0, 12)))
col1 = plot_grid(wunifa + theme(legend.position = "none"),
nrow = 1, ncol = 1, labels = c('A'),
label_size = 20)
col2 = plot_grid(wunioa + theme(legend.position = "none"),
nrow = 1, ncol = 1, labels = c('B'),
label_size = 20)
plot = plot_grid(col1, col2, legend1, nrow = 1, ncol = 3,
rel_widths = c(1.5, 1.5, 0.75), align = "h", axis = "t")
plot
dev.off()## quartz_off_screen
## 2Unweighted UniFrac
Fecal
adonis2(uuni_fecala ~ ethnicity, data = relab_fecal_filtered_after,
permutations = 4999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 4999
##
## adonis2(formula = uuni_fecala ~ ethnicity, data = relab_fecal_filtered_after, permutations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## ethnicity 1 0.20279 0.07361 2.3043 0.0044 **
## Residual 29 2.55215 0.92639
## Total 30 2.75494 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Oral
adonis2(uuni_orala ~ ethnicity, data = relab_oral_after,
permutations = 4999)## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 4999
##
## adonis2(formula = uuni_orala ~ ethnicity, data = relab_oral_after, permutations = 4999)
## Df SumOfSqs R2 F Pr(>F)
## ethnicity 1 0.04386 0.06782 1.0914 0.3346
## Residual 15 0.60285 0.93218
## Total 16 0.64671 1.00000NMDS visualization of unweighted UniFrac distances with after samples.
## Run 0 stress 0.1657305
## Run 1 stress 0.2285569
## Run 2 stress 0.1658235
## ... Procrustes: rmse 0.005360905 max resid 0.02528318
## Run 3 stress 0.1789747
## Run 4 stress 0.165719
## ... New best solution
## ... Procrustes: rmse 0.001621189 max resid 0.005848108
## ... Similar to previous best
## Run 5 stress 0.2010731
## Run 6 stress 0.1678215
## Run 7 stress 0.1824019
## Run 8 stress 0.182164
## Run 9 stress 0.1681404
## Run 10 stress 0.2237202
## Run 11 stress 0.1675407
## Run 12 stress 0.1678212
## Run 13 stress 0.223468
## Run 14 stress 0.1681404
## Run 15 stress 0.1781267
## Run 16 stress 0.1658233
## ... Procrustes: rmse 0.005827622 max resid 0.02552845
## Run 17 stress 0.2010731
## Run 18 stress 0.1657301
## ... Procrustes: rmse 0.001529493 max resid 0.005840392
## ... Similar to previous best
## Run 19 stress 0.1904371
## Run 20 stress 0.1657305
## ... Procrustes: rmse 0.001633068 max resid 0.005843246
## ... Similar to previous best
## *** Solution reached
## Run 0 stress 0.1364074
## Run 1 stress 0.1543538
## Run 2 stress 0.1367323
## ... Procrustes: rmse 0.04383167 max resid 0.1382546
## Run 3 stress 0.1304384
## ... New best solution
## ... Procrustes: rmse 0.04509471 max resid 0.1268594
## Run 4 stress 0.1384731
## Run 5 stress 0.1364076
## Run 6 stress 0.1385508
## Run 7 stress 0.1438585
## Run 8 stress 0.1543538
## Run 9 stress 0.1304384
## ... Procrustes: rmse 2.125227e-05 max resid 6.013072e-05
## ... Similar to previous best
## Run 10 stress 0.1304388
## ... Procrustes: rmse 0.0009354498 max resid 0.002662176
## ... Similar to previous best
## Run 11 stress 0.1364078
## Run 12 stress 0.1364078
## Run 13 stress 0.3689216
## Run 14 stress 0.1363188
## Run 15 stress 0.1304383
## ... New best solution
## ... Procrustes: rmse 0.000102872 max resid 0.0002883375
## ... Similar to previous best
## Run 16 stress 0.1694527
## Run 17 stress 0.1384728
## Run 18 stress 0.1304383
## ... New best solution
## ... Procrustes: rmse 0.0003611925 max resid 0.001027821
## ... Similar to previous best
## Run 19 stress 0.1380133
## Run 20 stress 0.1367322
## *** Solution reached
tiff(file = "figures/nmds_plot_taxa_uuni_after.tif", res = 150, width = 18,
height = 7, units = "in")
legend1 = get_legend(uunifa + theme(legend.box.margin = margin(0, 0, 0, 12)))
col1 = plot_grid(uunifa + theme(legend.position = "none"),
nrow = 1, ncol = 1, labels = c('A'),
label_size = 20)
col2 = plot_grid(uunioa + theme(legend.position = "none"),
nrow = 1, ncol = 1, labels = c('B'),
label_size = 20)
plot = plot_grid(col1, col2, legend1, nrow = 1, ncol = 3,
rel_widths = c(1.5, 1.5, 0.75), align = "h", axis = "t")
plot
dev.off()## quartz_off_screen
## 2Alpha Diversity
Alpha diversity was calculated using species-level metaphlan results in vegan (2.5-7). Shannon diversity and Shannon-Weaver evenness (Pielou’s evenness) were calculated.
shannond_fecal = diversity(relab_fecal_filtered[,2:496], index = "shannon")
shannond_oral = diversity(relab_oral[,2:496], index = "shannon")
shannone_fecal = shannond_fecal/log(specnumber(relab_fecal_filtered[,2:496]))
shannone_oral = shannond_oral/log(specnumber(relab_oral[,2:496]))
fecal_diversity = cbind(relab_fecal_filtered[,1], relab_fecal_filtered[,497:520], shannond_fecal, shannone_fecal)
fecal_diversity$Study_period = as.factor(fecal_diversity$Study_period)
oral_diversity = cbind(relab_oral[,1], relab_oral[,497:519], shannond_oral, shannone_oral)
oral_diversity$Study_period = as.factor(oral_diversity$Study_period)A linear model was used to examine differences between self-reported ethnicities and study period (before/during/after) in diversity.
div_fecal = lm(shannond_fecal ~ ethnicity*Study_period,
data = fecal_diversity)
summary(div_fecal)##
## Call:
## lm(formula = shannond_fecal ~ ethnicity * Study_period, data = fecal_diversity)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.3741 -0.2227 0.1292 0.3666 0.6455
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.825245 0.117163 24.114 <2e-16 ***
## ethnicityCaucasian 0.023833 0.180926 0.132 0.895
## Study_periodBefore -0.001497 0.190236 -0.008 0.994
## Study_periodDuring 0.117951 0.143495 0.822 0.413
## ethnicityCaucasian:Study_periodBefore -0.202191 0.293399 -0.689 0.492
## ethnicityCaucasian:Study_periodDuring 0.002998 0.222444 0.013 0.989
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4971 on 105 degrees of freedom
## Multiple R-squared: 0.03296, Adjusted R-squared: -0.01309
## F-statistic: 0.7157 on 5 and 105 DF, p-value: 0.613Anova(div_fecal)## Anova Table (Type II tests)
##
## Response: shannond_fecal
## Sum Sq Df F value Pr(>F)
## ethnicity 0.0023 1 0.0094 0.9229
## Study_period 0.7211 2 1.4591 0.2371
## ethnicity:Study_period 0.1599 2 0.3236 0.7243
## Residuals 25.9444 105summary(glht(div_fecal,linfct=mcp(Study_period="Tukey")))## Warning in mcp2matrix(model, linfct = linfct): covariate interactions found --
## default contrast might be inappropriate##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = shannond_fecal ~ ethnicity * Study_period, data = fecal_diversity)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## Before - After == 0 -0.001497 0.190236 -0.008 1.000
## During - After == 0 0.117951 0.143495 0.822 0.687
## During - Before == 0 0.119449 0.171249 0.698 0.763
## (Adjusted p values reported -- single-step method)even_fecal = lm(shannone_fecal ~ ethnicity*Study_period,
data = fecal_diversity)
summary(even_fecal)##
## Call:
## lm(formula = shannone_fecal ~ ethnicity * Study_period, data = fecal_diversity)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.31337 -0.03586 0.02678 0.07622 0.18643
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.649303 0.024716 26.271 <2e-16 ***
## ethnicityCaucasian 0.005016 0.038167 0.131 0.896
## Study_periodBefore 0.030051 0.040131 0.749 0.456
## Study_periodDuring 0.019407 0.030270 0.641 0.523
## ethnicityCaucasian:Study_periodBefore -0.086399 0.061893 -1.396 0.166
## ethnicityCaucasian:Study_periodDuring -0.003113 0.046925 -0.066 0.947
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1049 on 105 degrees of freedom
## Multiple R-squared: 0.03584, Adjusted R-squared: -0.01007
## F-statistic: 0.7807 on 5 and 105 DF, p-value: 0.5658Anova(even_fecal)## Anova Table (Type II tests)
##
## Response: shannone_fecal
## Sum Sq Df F value Pr(>F)
## ethnicity 0.00359 1 0.3266 0.5689
## Study_period 0.01186 2 0.5395 0.5846
## ethnicity:Study_period 0.02733 2 1.2427 0.2928
## Residuals 1.15454 105summary(glht(even_fecal,linfct=mcp(Study_period="Tukey")))## Warning in mcp2matrix(model, linfct = linfct): covariate interactions found --
## default contrast might be inappropriate##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: lm(formula = shannone_fecal ~ ethnicity * Study_period, data = fecal_diversity)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## Before - After == 0 0.03005 0.04013 0.749 0.732
## During - After == 0 0.01941 0.03027 0.641 0.796
## During - Before == 0 -0.01064 0.03613 -0.295 0.953
## (Adjusted p values reported -- single-step method)div_oral = lm(shannond_oral ~ ethnicity*Study_period, data = oral_diversity)
summary(div_oral)##
## Call:
## lm(formula = shannond_oral ~ ethnicity * Study_period, data = oral_diversity)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.32262 -0.10843 -0.04301 0.11353 0.66514
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.37870 0.06839 49.403 <2e-16 ***
## ethnicityCaucasian -0.05358 0.10658 -0.503 0.619
## Study_periodBefore -0.04956 0.09672 -0.512 0.612
## ethnicityCaucasian:Study_periodBefore 0.06253 0.15073 0.415 0.681
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2163 on 30 degrees of freedom
## Multiple R-squared: 0.01195, Adjusted R-squared: -0.08685
## F-statistic: 0.121 on 3 and 30 DF, p-value: 0.9471Anova(div_oral)## Anova Table (Type II tests)
##
## Response: shannond_oral
## Sum Sq Df F value Pr(>F)
## ethnicity 0.00410 1 0.0877 0.7691
## Study_period 0.00482 1 0.1030 0.7504
## ethnicity:Study_period 0.00805 1 0.1721 0.6812
## Residuals 1.40317 30even_oral = lm(shannone_oral ~ ethnicity*Study_period, data = oral_diversity)
summary(even_oral)##
## Call:
## lm(formula = shannone_oral ~ ethnicity * Study_period, data = oral_diversity)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.062819 -0.016022 -0.001432 0.020062 0.057567
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.724772 0.009848 73.599 <2e-16 ***
## ethnicityCaucasian -0.030771 0.015346 -2.005 0.054 .
## Study_periodBefore -0.008127 0.013926 -0.584 0.564
## ethnicityCaucasian:Study_periodBefore 0.031961 0.021703 1.473 0.151
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03114 on 30 degrees of freedom
## Multiple R-squared: 0.124, Adjusted R-squared: 0.03645
## F-statistic: 1.416 on 3 and 30 DF, p-value: 0.2574Anova(even_oral)## Anova Table (Type II tests)
##
## Response: shannone_oral
## Sum Sq Df F value Pr(>F)
## ethnicity 0.0018014 1 1.8576 0.1830
## Study_period 0.0002153 1 0.2221 0.6409
## ethnicity:Study_period 0.0021031 1 2.1687 0.1513
## Residuals 0.0290921 30Shannon diversity and evenness did not vary across self-reported ethnicities or study periods in either fecal or oral samples.
## quartz_off_screen
## 2## quartz_off_screen
## 2